function [t, h] = nonlin_sim

addpath ../

p = params;
sp = sim_par;

[x, z] = meshgrid(sp.x, sp.z);

c0 = 0.3 * ones(size(x));
h0 = 1;

X0 = [reshape(c0, sp.N * sp.Nz, 1); h0];

[t, X] = ode23s(@(t, X) pde(t, X, p, sp), [0 1e4], X0);


for i = 1:length(t)
    
    h(i) = X(i,end);
%     c = reshape(X(i,1:end-1), sp.Nz, sp.N);

    
    
end

ha = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);
semilogx(t, h, t, ha, 'k*');



function G = pde(t, X, p, sp)

h = X(end);
c = reshape(X(1:end-1), sp.Nz, sp.N);

F = zeros(size(c));
G = zeros(size(X));

eta = sp.z';

c_top = c(sp.Nz,:);
c_np1 = c(sp.Nz-1,:) - 2 * sp.dz * p.delta * h * c_top .* (1 - c_top);

% [u, w] = comp_velocity(c_top, h, p, sp);

u = zeros(sp.Nz, sp.N);
w = u;

dh = w(sp.Nz, 1) - p.delta * c(sp.Nz, 1);


% c_x = zeros(size(u));
% c_xx = c_x;

c_hat = fft(c')';
c_x = real(ifft(1i * sp.K .* c_hat, [], 2));
c_xx = real(ifft(-sp.K2 .* c_hat, [], 2));


% for j = 1:sp.Nz
%     
%     c_hat = fft(c(j,:));
%     c_x(j,:) = real(ifft(1i * sp.k .* c_hat));
%     c_xx(j,:) = real(ifft( -sp.k.^2 .* c_hat));
%     
% end

F(1,:) = c_xx(1,:) + 1 / h^2 * (2 * c(2,:) -2*c(1,:)) / sp.dz / sp.dz;


F(2:end-1,:) = (dh * eta(2:end-1) * ones(1,sp.N) - w(2:end-1,:)) / h .* (c(3:end,:) - c(1:end-2,:)) / 2 / sp.dz ...
    - u(2:end-1,:) .* c_x(2:end-1,:) + c_xx(2:end-1,:) + ...
    (c(3:end,:) - 2 * c(2:end-1,:) + c(1:end-2,:)) / sp.dz^2 / h^2;


F(end,:) = (p.delta * c_top).^2 .* (1 - c_top) - u(end,:) .* c_x(end,:) + c_xx(end,:) ...
    + (c_np1 - 2 * c(end,:) + c(end-1,:)) / sp.dz^2 / h^2;
        

G(1:end-1) = reshape(F, sp.N * sp.Nz, 1);
G(end) = dh;